###########################################################
##
## Distance to Lake Chad
##
###########################################################

today <- Sys.Date()

source(paste0(wkdir,'/udf/regional_seas.R'))

##
## Input
##

## fishnet X adm0 boundary N=33252
fishnet_lake_chad_ntl_adm0_x_d01_shp <- paste0(wkdir,'/local/gis_data/003_boundaries/fishnet/fishnet_lake_chad_ntl_adm0_x_d01')
fishnet_adm0_only <- st_read(dirname(fishnet_lake_chad_ntl_adm0_x_d01_shp),basename(fishnet_lake_chad_ntl_adm0_x_d01_shp))
st_crs(fishnet_adm0_only) <- "EPSG:4326"

fishnet_adm0_only.prj <- fishnet_adm0_only %>%
  # 102022
  st_transform(crs="+proj=aea +lat_1=20 +lat_2=-23 +lat_0=0 +lon_0=25 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m no_defs") %>% 
  st_centroid() %>% 
  st_geometry()

M.fishnet.zones <- fishnet_adm0_only.prj %>% 
  st_transform(4326) %>%
  st_coordinates() 

# convert polygons to points 
fishnet.zones.pts = st_as_sf(as.data.frame(M.fishnet.zones), coords = c("X", "Y"), crs = 4326)
fishnet.zones = st_drop_geometry(subset(fishnet_adm0_only,select=c(OBJECTID, GID_0)))
fishnet.zones.pts = dplyr::bind_cols(fishnet.zones.pts,fishnet.zones)

##
## LAKE CHAD BOUNDARY
## OLD
##

fishnet.zones.pts$dist2lakechadold<-NA
big_lake <- st_read(paste0(wkdir,'/local/gis_data/Lakes_for_Brian/big_lake.shp'))
st_crs(big_lake)
fishnet.zones.pts$dist2lakechadold <- unlist(nngeo::st_nn(fishnet.zones.pts, big_lake, k = 1, returnDist = T)$dist) / 1000


##
## LAKE CHAD BOUNDARY
## NEW
##
fishnet.zones.pts$dist2lakechadnew<-NA
ne_10m_lakes <- st_read(paste0(wkdir,'/local/gis_data/Lakes_for_Brian/ne_10m_lakes.shp'))
st_crs(ne_10m_lakes)
fishnet.zones.pts$dist2lakechadnew <- unlist(nngeo::st_nn(fishnet.zones.pts, ne_10m_lakes, k = 1, returnDist = T)$dist) / 1000
summary(fishnet.zones.pts$dist2lakechadnew)

dist2lakechad_fn = file.path(wkdir,"proc_data","DIST2LAKECHAD.dta")
dist2lakechad_out <- fishnet.zones.pts %>%
  sf::st_drop_geometry() %>%
  dplyr::select(OBJECTID,GID_0,dist2lakechadold,dist2lakechadnew) %>%
  dplyr::rename_with(.,tolower) 
names(dist2lakechad_out)

# Adding the label value
attr(dist2lakechad_out$dist2lakechadold, "label") <- "Euclidean distance to former large Lake Chad : Source?"
attr(dist2lakechad_out$dist2lakechadold, "label")
attr(dist2lakechad_out$dist2lakechadnew, "label") <- "Euclidean distance to current (small) Lake Chad : Natural Earth"
attr(dist2lakechad_out$dist2lakechadnew, "label")
attr(dist2lakechad_out, "label") <- sprintf("dist2lakechad.R last run on %s",Sys.Date())
head(dist2lakechad_out)
summary(dist2lakechad_out)

## Save the data
haven::write_dta(dist2lakechad_out, dist2lakechad_fn)